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Two  problems  were  addressed  during  this  effort:  numerical  diffusion, 

which  exists  in  HULL,  significantly  affects  the  calculations  of  nuclear  air 
bursts  a  low  overpressure,  and  (Zf  some  HULL  calculations  being  performed  are 
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ity  in  them  and  to  efficiently  perform  these  types  of  calculations,  modifica¬ 
tions  to  HUE!,  were  needed. 

Modifications  were  made  to  the  HULL  code  in  an  attempt  to  reduce  numerical 
diffusion  at  low  overpressure  by  employing  a  version  of  Flux  Corrected  . - 
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v Transport  (FCT).  A  version  of  FCT  was  provided  by  the  AFWL  and  this  was 
implemented  in  HULL  and  tested  for  a  simulation  of  a  spherically  symmetric 
explosion.  The  modifications  to  HULL  were  made  in  the  AFWL  two-dimensional 
version. 

The  low  overpressure  shocks  were  sharpened  by  the  FCT  approach  for  a  given 
cell  size.  The  improvement  is  not  considered  to  be  sufficient  to  justify  the 
additional  cost  incurred  as  a  result  of  the  more  complex  numerical  differenc¬ 
ing  in  this  particular  version  of  FCT. 

In  the  second  area,  a  technique  was  developed  for  avoiding  computing  within 
certain  regions  of  the  mesh,  which  was  also  implemented  in  the  two-dimensional 
version,  and  additional  quantitative  criteria  were  established  as  the  basis 
for  subsequent  development  of  better  rezone  techniques.* 

An  implementation  of  a  computational  subgrid  was  developed  and  is  recom¬ 
mended  for  use  in  HULL  to  obviate  the  need  for  rezone  fn\many  problems.  In 
addition,  it  was  found  that  numerical  errors  in  HULL  associated  with  zone  size 
changes  are  less  pronounced  than  the  effects  introduced  with  different 
artificial  viscosity  prescriptions.  Furthermore,  it  was  observed  that  pulse 
dispersion  introduced  by  propagation  through  regions  of  coarse  zoning  will 
lead  to  reductions  in  peak  overpressure,  which  in  regions  of  low  overpressure 
will  never  be  regained  as  is  the  case  at  high  overpressure. 
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The  research  reported  herein  was  directed  towards  making  some 
improvements  to  the  HULL  code.  An  AFWL  version  of  the  Flux  Corrected 
Transport  (FCT)  was  implemented  in  HULL  and  compared  with  the  usual 
HULL  difference  technique  for  an  explosion  of  a  high  pressure  isothermal 
sphere  into  air.  In  addition,  algorithms  for  treating  only  a  portion 
of  the  computational  mesh,  a  computational  subgrid,  were  developed. 

One  dimensional  HULL  calculations  were  performed  for  use  in  the  develop¬ 
ment  of  improved  rezone  techniques. 
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the  computation  subgrid  technique.  Dr.  Willard  R.  Thomas  supervised 
the  one-dimensional  HULL  calculations.  The  authors  especially  thank 
Ms.  Caroline  Peerson-Reeves ,  who  provided  technical  assistance  to 
Drs.  Hasdai  and  Thomas. 
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SECTION  1 


INTRODUCTION 

The  effort  reported  herein  was  performed  to  complement  the 
Defense  Nuclear  Agency  (DNA)/Air  Force  Weapons  Laboratory  (AFWL)  pro¬ 
gram  addressing  aircraft  survivability  during  base  escape.  The  AFWL 
was  sponsored  by  DNA  to  perform  fine-zoned  HULL  calculations  in  order 
to  provide  an  improved  representation  of  air  blast  at  ground  level  for 
low  overpressures  (several  psi)  as  well  as  air  blast  environments  at 
different  heights-of-target.  These  data  are  required  for  higher  con¬ 
fidence  assessments  of  aircraft  survi vabi 1 i ty  while  aircraft  (bombers 
and/or  tankers)  are  escaping  the  effects  of  a  hostile  nuclear  attack. 

The  principal  requirement  for  the  AFWL  calculations  is  to  reduce  the 
uncertainties  in  air  blast  at  low  overpressures .  Although  much  work 
has  been  done  in  this  area  considerable  uncertainty  still  exists.  In 
particular,  the  uncertainties  in  range  to  which  certain  overpressures 
extend  are  making  base  escape  assessments  difficult. 

SAI  conducted  two  tasks  in  support  of  the  above  program.  The 
first  was  to  implement  an  AFWL  version  of  the  Flux  Corrected  Transport 
(FCT)  technique  developed  by  Boris,  et  al  (Reference  1-3)  into  HULL  to 
allow  better  treatment  of  low  overpressure  shocks.  The  second  was  to 
develop  improved  rezone  techniques  to  be  used  in  HULL.  A  rezone  capa¬ 
bility  is  used  to  maintain  adequate  resolution  of  shocks  without  requir¬ 
ing  many  zones  in  regions  where  little  hydrodynamic  activity  is  occur¬ 
ring. 
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SECTION  2 


THE  HULL  HYDROCODE 

HULL  is  the  acronym  for  a  computational  hydrodynamics  code 
used  to  solve  the  nonlinear  hyperbolic  conservation  laws  for  an  ideal 
inviscid  fluid.  It  exists  in  many  versions  (references  4,5),  any  one 
of  which  may  be  invoked  when  a  calculation  is  to  be  performed.  The 
three-dimensional  version  of  HULL,  for  example,  solves  finite  dif¬ 
ference  representations  of  a  closed  system  of  partial  differential 
equations.  These  equations  describe  nondisippative  continuum  fluid 
flow  in  the  absence  of  electric  and  magnetic  fields.  The  equations 
include  the  conservation  law  for  mass,  momentum,  and  energy,  written 
in  a  Lagrangian  frame,  as  well  as  the  fluid  equation  of  state,  and 
are  listed  here  for  direct  reference. 


do 

dt 

pV-TT  =  0 

(D 

diT 

pdt  + 

VP  =  pg" 

(2) 

dE  , 
pdt  + 

V-Pu  =  puT 

(3) 

P  = 

P  (p,I) 

(4) 

where  .  is  material  density,  P  is  pressure,  u  is  fluid  velocity,  E  is 
total  specific  energy,  I  is  internal  specific  energy,  g"  is  acceleration 
due  to  gravity,  and  t  is  time.  Although  the  equations  are  written  in  a 
Lagrangian  frame,  during  one  phase  of  the  hydrocode  calculation,  masses. 
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momentum,  and  energy  are  transported  across  a  fixed  calculational  grid 
making  HULL  an  Eulerian  code.  A  first  order  donor  cell  technique  is 
used;  this  introduces  artificial  increases  in  entropy  that  usually  per¬ 
mit  the  calculation  of  strong  shock  waves  without  the  introduction  of 
a  physical  representation  for  fluid  viscosity. 

Although  additional  information  on  HULL  are  available  in  the 
references  cited  above.  Appendix  A  is  provided  to  acquaint  the  reader 
with  the  capability  of  HULL  to  accurately  predict  airblast  from  HE 
detonations.  HULL  has  been  demonstrated  to  be  very  accurate  for  many 

9 

HE  events,  however  considerable  uncertainty  still  exists  for  the  low 
overpressure  regime.  For  example,  in  the  Appendix  the  discussion  of 
Figure  27  (comparisons  of  overpressure  measurements  on  Dial  Pack  with 
SHELL  results)  points  out  that  results  from  SHELL  below  10  psi  fall 
below  the  data.  The  use  of  HULL  improved  this  result  somewhat;  however, 
there  exists  considerable  room  for  improvement  for  the  several  psi 
regime.  From  the  figures  in  Appendix  A  it  is  clear  that  there  can 
exist  differences  between  40  and  807  in  range  to  a  given  low  overpres¬ 
sure. 

This  large  uncertainty  in  range  leads  to  significant  differences 
in  aircraft  system  survi vabi 1 i ty  assessments.  This  effort  addressed 
two  areas  where  Chandler  and  Needham  of  AFWL  felt  improvement  in  HULL 
was  desirable  and  would  help  reduce  these  large  uncertainties.  Chandler 
felt  improvements  could  be  gained  by  including  FCT  methods,  originally 
developed  at  NRL  (references  1-3),  in  the  HULL  code.  He  performed 
some  work  in  one-dimension  (reference  6)  which  demonstrated  that  FCT 
techniques  can  improve  shock  definition  when  used  with  the  HULL  dif¬ 
ference  technique.  SAI  implemented  the  AFWL  FCT  technique  in  HULL  in 
two  dimensions  and  then  tested  it. 

The  second  task  was  performed  because  of  a  need  for  improving 
the  rezone  techniques  identified  by  Needham.  Rezone  techniques,  for 
the  purposes  of  this  effort,  attempt  to  provide  fine  zoning  where 


shocks  exist  and  coarse  zoning  where  little  hydrodynamic  activity  is 
occurring.  The  intent  is  to  reduce  the  number  of  zones  used  in  the 
calculations  to  an  optimum  value,  thereby  obtaining  the  most  informa¬ 
tion  for  the  least  cost.  However,  if  used  improperly,  a  rezone  technique 
may  reduce  the  number  of  zones  in  regions  that  in  fact  have  significant 
hydrodynamic  activity. 
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SECTION  3 


IMPLEMENTATION  OF  FCT  IN  HULL 

To  support  implementation  of  a  flux  corrected  transport  (FCT) 
method,  SAI  obtained  the  FCT  difference  equations  from  AFWL,  designed 
the  architectural  modifications  to  HULL,  wrote  the  software  to  imple¬ 
ment  these  modifications,  tested  the  software  in  a  simplified  version 
of  HULL,  and  then  provided  the  coding  to  AFWL. 

The  small,  self-contained  version  of  HULL  developed  by  SAI  for 
timing  studies  on  the  CRAY-1  computer  (reference  7)  was  chosen  as  the 
version  for  the  initial  implementation  of  FCT  for  several  reasons,  the 
more  important  of  which  are:  (1)  the  code  architecture  is  similar  to 
production  versions  of  HULL,  thus  the  FCT  code  will  not  be  unique  to 
this  version,  (2)  execution  is  significantly  faster,  thereby  minimizing 
cost  and  turnaround;  also  most  of  the  HULL  system  bookkeeping  overhead 
is  absent  for  this  in-core  version,  (3)  this  version  will  execute  on 
the  CDC  machines  available  at  AFWL.  The  method  of  FCT  implementation 
employed  is  presented  following  a  short  discussion  of  the  computational 
phases  concept  used  in  HULL. 

3.1  HULL  PHASES 

The  fundamental  variables  computed  and  stored  in  HULL  two- 
dimensional  airblast  calculations  are  pressure,  two  components  of 
velocity,  internal  energy,  and  mass  of  each  cell  of  the  mesh  of  cells 
which  represent  the  physical  conditions  simulated  by  a  calculation. 

In  many  calculations,  the  total  of  the  mesh  variables  exceeds  central 
memory  allocations  of  most  machines  for  which  HULL  is  implemented, 
e.g.,  a  100  by  200  cell  single  material  mesh  requires  100,000  storage 
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locations  just  for  the  mesh  variables.  Therefore,  mesh  variables 
must  generally  be  placed  in  auxiliary  storage,  currently  on-line  disk 
or  extended  core.  The  mesh  variables  are  stored  in  these  media  in  an 
array  which  corresponds  to  a  geometric  image  of  the  physical  conditions. 
In  two  dimensional  calculations,  cell  indexes  (I,J)  correspond  to 
geometric  coordinates  (X,Y).  A  row  of  a  two  dimensional  mesh  is  that 
set  of  variables  for  which  J  =  constant  (Y  =  constant)  and  1  <  I  <  IMAX 
(XMIN  5  X  ^  XMAX).  In  this  storage  mode,  contiguous  rows  of  the  mesh 
are  sequentially  accessible. 

The  computational  phases  in  HULL  are  distinct,  independent 
numerical  operations.  A  complete  set  of  these  phases  operating  in 
sequence  will  advance  the  mesh  from  time  t  to  time  t  +  dt.  This  is 
called  a  time  step.  Repetitive  operation  of  this  set  upon  the  mesh 
produces  the  time  evolution  of  the  hydrodynamic  variables  (mesh)  in  a 
calculation.  The  computational  phases  are  independent  in  the  sense 
that  no  intermediate  communication  occurs  between  the  phases  other 
than  by  means  of  the  mesh  variables. 

Combining  these  storage  and  phasing  methods,  a  successful 
code  structure  would  be  that  of  operating  upon  the  entire  mesh  in 
turn  with  each  phase.  This  structure  requires  complete  input/output 
processing  of  the  mesh  for  each  phase,  unless  the  mesh  can  be  contained 
in  central  memory.  The  actual  numerical  computations  used  in  the  HULL 
differencing  require  at  most  two  rows  for  a  phase  plus  space  for  tem¬ 
porary  storage  which  is  internal  to  the  phase.  Figure  1  illustrates 
the  phase  set  operation  structure  used  in  HULL  which,  in  contrast  to 
the  above  mentioned  structure,  requires  one  input/output  processing 
of  the  mesh  per  time  step.  The  description  of  the  phases  indicated  in 
Figure  1  is: 
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HULL  phases  and  phase  operation  sequence. 


•  EOS  -  Equation  of  State,  updates  cell  pressure  to 

time  t  value. 

•  HYDRO  -  Lagrangian  phase,  computes  intermediate 

internal  energy  and  velocity  component 
values  on  a  Lagrangian  mesh  at  time  t  +  dt. 

•  FLUX  -  Transport  phase,  fluxes  mass  and  internal 

energy  between  cells  to  obtain  final  mesh 
variable  values  at  time  t  +  dt. 

The  separate  phases  are  provided  a  pointer  index  indicating 
the  row  to  be  processed  by  each.  The  phases  are  executed  in  the  order 
EOS,  HYDRO,  FLUX  -  which  has  the  effect  that  each  phase  in  i ffect 
requires  one  row  of  the  mesh  in  core;  i.e.,  that  row  updated  by  that 
phase.  In  particular,  in  terms  of  Figure  1,  after  the  EOS  operation 
upon  row  J,  the  HYDRO  phase  can  update  row  J-l  variables  using  mesh 
values  at  time  t  from  rows  J  and  J-l.  Upon  completion  of  HYDRO,  FLUX 
can  produce  the  time  t  +  dt  value  of  the  variables  in  row  J-2  using 
the  intermediate  mesh  values  from  rows  J-l  and  J-2.  Rows  J+l  to  JMAX, 
and  rows  1  to  J-3  are  not  involved  in  this  phase  set  operation  and 
need  not  be  in  central  memory.  Prior  to  the  next  phase  set  operation 
the  pointer  indexes  are  advanced  one  row,  row  J-2  now  at  time  t  +  dt 
is  output,  and  row  J  +  1  at  time  t  is  input.  After  each  row  has  been 
operated  upon  by  each  phase,  the  entire  mesh  has  been  advanced  to  time 
t  +  dt. 

It  is  evident  that  this  structure  is  independent  of  the  particu 
lar  storage  used  for  those  rows  of  the  mesh  not  currently  designated  by 
the  phase  row  pointer  indexes.  It  is  compatible  with  central  memory, 
extended  core,  or  disk  storage  of  the  mesh. 
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3.2  FLUX  CORRECTED  TRANSPORT  PHASING 

The  Flux  Corrected  Transport  (FCT)  algorithm  as  it  is  currently 
configured  by  AFWL ,  can  be  separated  into  two  independent  numerical 
operations,  termed  diffusion  (DIFF)  and  anti-diffusion  (ADIFF).  Further, 
these  operations  can  be  designated  as  operational  phases,  which  are 
quite  similar  in  structure  to  the  HULL  phases.  This  approach  is  the 
one  which  has  been  taken  for  implementation  of  FCT  in  HULL,  in  which 
the  advantages  of  the  HULL  phasing  structure  are  apparent  in  minimiza¬ 
tion  of  storage  requirements. 

The  FCT  algorithm  operates  upon  the  mesh  variables  in  conserved 
quantity  form;  i.e.,  mass,  individual  momentum  components,  and  total 
energy.  The  available  algorithm  is  constructed  for  two  spatial  dimen¬ 
sion  calculations,  which  is  the  only  form  implemented  during  this 
effort.  There  are  then,  four  conserved  quantities  per  cell  to  consider. 

The  FCT  phases  are  independent  of  each  other  in  the  sense  that 
the  HULL  phases  of  Figure  1  are  independent.  However,  the  FCT  phases 
are  not  independent  of  the  HULL  phases  in  this  sense.  They  are  not 
independent,  because  each  one  requires  information  other  than  that  con¬ 
tained  in  the  mesh  variables.  The  effects  of  non-independence  are  min¬ 
imized,  in  terms  of  required  extra  variable  storage,  by  sequencing  the 
combined  HULL  and  FCT  phase  set  as  is  illustrated  in  Figure  2.  When 
ADIFF  is  completed  upon  row  J-5,  this  row  can  be  output;  it  is  not 
needed  for  the  ADIFF  operation  upon  row  J-4  to  be  done  in  the  next 
operation  of  the  phase  set.  Therefore,  FCT  when  included  in  HULL  re¬ 
quires  six  rows  of  mesh  in  central  memory.  The  extra  storage  require¬ 
ments  arising  from  the  non- independence  are  defined  in  the  following 
discussion  of  the  FCT  phases. 
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HULL  and  FCT  phases,  and 
phase  operation  sequence. 


3.3 


FCT  DIFFUSION  PHASE  -  DATA  STORAGE  REQUIREMENTS 

The  diffusion  of  the  conserved  quantity  variables  of  a  cell 
requires  the  following: 

•  The  conserved  quantities  of  the  subject  cell  and  those 
of  its  four  adjacent  cells  at  time  t, 

•  The  volumes  of  these  five  cells, 

•  The  boundary  velocities  on  the  four  boundaries  of  the 
subject  cel  1 , 

•  The  conserved  quantities  of  the  subject  cell  at  time 
t  +  dt. 

Figure  3  illustrates  the  geometric  arrangement  of  the  cells  of  the  mesh 
involved,  when  the  center  cell  (C)  is  undergoing  diffusion.  In  terms 
of  Figure  2,  cell  C  in  Figure  3  is  in  row  J-3  and  cells  A  and  B  are  in 
rows  J-2  and  J-4  respectively.  The  HYDRO  phase  indicated  in  Figure  2 
will  alter  the  mesh  variables  of  row  J-l  from  their  values  at  time  T. 
Therefore,  it  is  apparent  that  separate  storage  of  four  rows  of  the 
mesh  in  conserved  quantity  form  is  required  containing  the  values  of 
rows  J-l  through  J-4  at  time  T.  This  will  satisfy  the  first  data 
requirement  listed  above.  It  is  to  be  noted  that  the  EOS  phase  affects 
only  the  pressure  variables  in  the  mesh,  which  are  not  used  in  FCT;  and 
row  J  is  not  required  in  separate  storage. 

The  volume  of  each  cell  can  be  computed  from  the  cell  geometry 
arrays  DX ( I ) ,  DY(J),  and  TAU(I). 

TAU(I)  contains  the  area  between  the  concentric  circles  bounding 
each  cell  in  cylindrical  geometry,  of  which  the  mesh  is  a  vertical  cross- 
section  of  a  right  cylinder.  These  arrays  are  extant  in  HULL  and  addi¬ 
tional  storage  is  not  required.  The  second  data  requirement  is  then 
sati sf ied . 
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Figure  3.  DIFF  phase  mesh  and  velocity  configuration 


The  third  data  requirement  -  boundary  velocities  -  will  re¬ 
quire  separate  storage.  These  velocities  are  computed  as  temporary 
variables  internal  to  the  FLUX  phase.  The  velocities  are  extracted 
as  each  is  computed  in  FLUX,  and  placed  in  storage  by  rows.  These 
same  velocities  are  used  in  the  anti-diffusion  phase  of  FCT  and  the 
total  storage  requirements  are  established  in  the  discussion  of  ADIFF. 

The  center  cell  conserved  quantity  variables  at  timet  +  dt 
can  be  computed  directly  from  the  mesh.  After  the  FLUX  phase,  the 
mesh  variables  (mass,  velocity  components,  and  internal  energy)  are 
the  time  t  +  dt  values.  Therefore,  no  additional  storage  is  necessary 
for  the  fourth  data  requirement  of  the  DIFF  phase. 

The  diffused  row  of  the  mesh  produced  by  DIFF  can  be  stored  in 
the  four  corresponding  variables,  cell  by  cell  in  row  J-3  of  the  mesh. 
However,  this  would  result  in  the  situation  wherein  the  ADIFF  phase  must 
access  its  input  variables  from  two  different  storage  structures  -  an 
awkward  situation  at  best.  Thus  the  DIFF  phase  results  are  placed  in 
a  temporary  storage  array,  which  contains  the  five  rows  of  diffused 
conserved  quantity  variables  required  in  ADIFF. 

3.4  FCT  DIFFUSION  PHASE  -  ALGORITHM 

Figure  4  contains  the  form  of  the  velocity  function  used  in 
the  diffusion  phase.  The  functions  in  this  exhibit  are  in  terms  of 
the  representative  storage  array  of  Figure  3,  in  which  the  velocity 
function  values  are  assigned  indexable  locations  between  their  related 
cells.  The  form  of  the  boundary  velocities  (i.e.,  E  in  Figure  4)  is 
the  same  in  DIFF  and  ADIFF;  therefore,  each  velocity  is  converted  to 
the  variable  form  E  as  it  is  extracted  to  storage  from  the  FLUX 
phase. 


In  terms  of  Figure  3,  if  C  e  (1,0)  and 

L  e  (1-2,0),  etc. 

UFL  =  U F ( I - 1  ,J)  UFR  =  UF ( 1+1  ,0) 

UFB  =  UF( I , J-I )  UFA  =  UF( 1 ,0+1 ) 

where 

UF(i,j)  =  (l/6)-(l/2)E(i ,  j )  +  (l/3)(E(i,j))2 
E(i,j)  =  |  DT-VELOCITY  (i , j  )/DR( i ,  j )  | 

DT  =  (T+dT)-T  (TIME  STEP  ADVANCE) 

DR  =  AVERAGE  DISTANCE  BETWEEN  CELL  CENTERS 
VELOCITY  =  BOUNDARY  VELOCITY  FROM  FLUX  PHASE 

THUS: 

E  ( i  - 1  ,  J )  =  |  2  •PJjV.ELOCIJ.Y  _(_( 1-2,0)  XM)i| 

DX ( 1-2,0)  +  DX ( 1 ,0 ) 

E  ( I  +  1  ,0 )  =  1  2  •-DT/ VELOCITY  ((l  ,0)  >_(_I +2,0  )JJ 
DX ( I  ,0 )  +  DX ( 1+2 ,0 ) 

E ( 1 ,0-1  )  =  I  2  -DT- VELOCITY  ( (1,0-2)  ►  (I ,0 )  )J 
DY( I ,0-2)  +  DY ( I ,0 ) 

E(I,0+1)  =  1 2-DT- VELOCITY  ((1,0)  .  (.1,0+2  ),)J 
DY ( 1 ,0 )  +  DY ( I ,0+2 ) 
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Figure  5  contains  the  computational  form  of  the  diffusion 
phase  of  the  FCT  algorithm.  The  quantity  DQC  of  the  figure  (dif¬ 
fused  conserved  quantity)  is  obtained  for  each  conserved  quantity  of 
every  cell  in  the  row  upon  which  DIFF  is  operating,  row  J-3  of 
Figure  2. 

The  boundary  conditions  for  diffusion  are  relatively  straight¬ 
forward.  The  guiding  principle  in  establishing  the  computational 
rules  at  the  boundaries  is  that  the  diffusion  flux  of  the  conserved 
quantity  is  zero  across  the  mesh  boundaries,  and  is  zero  at  the  limits 
of  the  (I,J)  index  range  which  is  being  processed  by  FCT.  The  particu¬ 
lar  computational  method  chosen  is  to  set  the  appropriate  velocity 
function  equal  to  zero.  A  zero  velocity  does  not  yield  this  result, 
as  is  apparent  from  the  form  of  UF  (i,j)  in  Figure  4. 

3.5  FCT  ANTI-DIFFUSION  PHASE  DATA  STORAGE  REQUIREMENTS 

The  anti-diffusion  of  the  diffused  quantity  variables  of  a  cell 
requires  the  following: 

•  The  diffused  conserved  quantities  of  the  subject  cell 
and  those  of  the  twelve  third-nearest  neighbor  cells, 

•  The  volumes  of  these  thirteen  cells, 

•  The  boundary  velocities  of  the  sixteen  common  cell 
boundaries  of  the  thirteen  cells. 

Figure  6  illustrates  the  geometric  arrangement  of  the  above  quantities, 
when  the  center  cell  (C)  is  undergoing  anti-diffusion.  In  terms  of 
Figure  2,  cell  C  in  Figure  6  is  in  row  J-5  and  cells  AA  and  BB  are  in 
rows  J-3  and  J-7  respectively.  Recall  that  rows  J-3  through  J-7  are 
rows  of  diffused  conserved  quantities  which  are  in  a  storage  area 
separate  from  the  mesh,  thus  ADIFF  has  already  processed  mesh  rows 
J-6  and  J-7  which  are  not  required  by  ADIFF  and  need  not  be  in  central 
memory.  The  five  rows  of  diffused  conserved  quantity  variables  satisfy 
the  first  data  requirement  listed  above. 


In  terms  of  Figure  3,  if  C  =  ( I , J )  and 

L  =  (1-2, J),  etc. 

DQC  =  QC  +  (GFR-GFL)  +  (GFA-GFB) 

where 

DQC  -  DIFFUSED  CONSERVED  QUANTITY  OF  CELL  (I.J) 

QC  -  CONSERVED  QUANTITY  OF  CELL  (I.J)  AT  TIME  T+DT 

GF*  -  DIFFUSION  FLUX  OF  QC  ACROSS  CELL  (I.J)  BOUNDARIES, 
WHERE  *  IS  L,  R,  A,  B. 

and 

GF*  =  GF  (UF*.  QC,  Q*) 

GFL  =  (VOLUME  (1-2,  J)  +  VOLUME  ( I , J ) )( 1/2)  ’UFL • ( RC-RL) 
GFR  =  (VOLUME  (I.J)  +  VOLUME  ( I+2.J) ) (1/2) -UFR- (RR-RC) 
GrA  =  (VOLUME  (I.J+2)  +  VOLUME  (I.J))  ( 1/2) -UFA • (RA-RC) 
GFB  =  (VOLUME  (I.J-2)  +  VOLUME  (I.J))  (1/2) -UFB • (RC-RB) 

where 

UF*  -  ARE  AS  SPECIFIED  IN  FIGURE  4 

R*  -  CONSERVED  QUANTITY  DENSITY  -  Q*/VOLUME(*) 

WHERE  *  IS  C.A.B.L.R;  AND  Q*  IS  THE 
CONSIRVED  QUANTITY  AT  TIME  T. 


Tiyure  5.  DIFF  phase  -  conserved  quantity 
diffusion  computation. 
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REFER  TO  FIGURE  3  FOR  NOTATION  CONVENTION:  NOTING  THAT  LL  -  LEFT  LEFT. 
VF**IS  THE  BOUNDARY  VELOCITY  FUNCTION,  *  is  L,A,R,B. 

Figure  6.  ADIFF  phase  mesh  and  velocity  configuration. 
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The  boundary  velocity  array  storage  convention  is  that  the 
velocities  found  with  index  set  (I,J)  for  example  for  cell  B,  are  in 
the  relation  to  B  as  are  VFBL  and  VFBB  of  Figure  6.  Therefore,  the 
number  of  rows  of  boundary  velocities  required  in  storage  is  six,  be¬ 
cause  the  Y-direction  velocities  obtained  from  FLUX  operating  upon 
row  J-2  of  Figure  2  are  associated  with  the  below  boundary  of  row  J-l . 
There  are  two  boundary  velocities  per  cell,  so  that  the  six  rows  of 
storage  are  equivalent  to  three  rows  of  conserved  quantity  storage. 

The  cell  volumes,  as  in  DIFF,  can  be  computed  from  the  existing  cell 
geometry  arrays  and  require  no  additional  storage. 

The  resultant  conserved  quantities  of  FCT  produced  at  the  end 
of  ADIFF,  on  a  cell  by  cell  basis  are  converted  to  the  standard  repre¬ 
sentations  of  the  hydrodynamic  variables  and  replace  the  time  t  +  dt 
values  of  mass,  velocity  components,  and  internal  energy  in  row  J-5. 
Thus,  no  additional  storage  is  required  to  hold  the  fully  corrected 
conserved  quantities.  Also,  if  the  calculation  is  mul ti -material ,  the 
individual  material  masses  can  be  updated  by  ADIFF  to  sum  to  the  new 
corrected  mass  of  each  cell. 

3.6  FCT  ANTI-DIFFUSION  -  ALGORITHM 

Figure  7  contains  the  form  of  the  velocity  function  used  in 
the  anti-diffusion  phase.  The  functions  in  this  exhibit  are  in  terms 
of  the  storage  array  of  Figure  6,  in  which  the  velocity  functions  are 
assigned  locations  capable  of  being  indexed.  The  storage  array  of 
Figure  6  was  incorporated  into  the  coding  of  the  ADIFF  phase  sub¬ 
routine.  As  will  be  seen  from  the  description  of  the  corrector  func¬ 
tions,  this  type  of  array  is  accessible  by  each  part  of  the  corrector 
functions  directly  under  indexing  control. 


In  terms  of  Figure  6,  if  C  =  (I,J)  and 

LL  =  (1-4, J),  etc. 

VFL  =  V  F ( I - 1 ,J)  VFRR  =  VF(I+3,J) 

VFA  -  VF(I,J+1)  VFBL  =  VF(I-1  ,J-2) 

where 

VF(i,j)  =(l/6)(l-(E(i,j))2) 

E ( i , j )  =  [ DT* VELOC ITY( i ,j )/DR(i ,j) | ,  c.f.  Figure  4. 

Thus: 

E ( I - 1 , J -2 )  =  12-PT-VELOCITY  ((I-2.J-2)  -  (I,J-2))[ 

DX( I-2.J-2)  +  DX(I,J-2) 

v t  iii  =  |2-DT. VELOCITY  ((I+2.J)  -  (I+2,J+2))| 

OY ( 1+2, J )  +  DY( 1+2, J+2) 

The  remaining  fourteen  values  of  E  are  constructed  in  a  similar  manner. 


Figure  7.  ADIFF  phase  velocity  function. 
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Figure  8  contains  the  computational  form  of  the  anti¬ 
diffusion  phase  of  the  FCT  algorithm.  The  quantity  FQC  (flux  cor¬ 
rected  conserved  quantity)  is  obtained  for  each  conserved  quantity 
of  every  cell  in  the  row  upon  which  ADIFF  is  operating,  row  J-5  of 
Figure  2.  The  formalism  of  DIFF  and  ADIFF  is  quite  similar  with  the 
exception  of  operation  of  the  corrector  function  (K  of  Figure  8) 
upon  the  diffused  conserved  quantity  fluxes  before  a  corrected  con¬ 
served  quantity  is  computed.  The  corrector  function  is  presented  in 
the  following  discussion. 

The  corrector  function  is  a  function  of  the  value  and  of  the 
sign  of  the  flux  upon  which  it  is  operating.  This  dependence  separates 
the  functional  form  into  three  routines  which  are  presented  in  the 
order  of  operation  employed  in  the  coding.  Figure  9  contains  the 
computational  form  of  the  part  of  the  corrector  function  termed  Test  1. 
This  test  employs  the  values  of  the  conserved  quantities  of  cells 
(1-4, J)  to  (1+4, J)  and  (I,J-4)  to  (I , J+  4 )  in  terms  of  Figure  6,  or 
the  in-line  third  and  first  nearest  neighbor  cells.  This  test  of 
itself  requires  that  five  rows  of  diffused  conserved  quantities  be 
available  to  ADIFF.  Step  1  of  Test  1  is  self-evident  from  the  form  of 
the  GF*  function  of  Figure  8  and  Figure  5  ;  it  is  the  result  of  the 
fact  that  the  adjacent  cells  have  identical  values  of  the  conserved 
quantity  being  corrected.  Parts  2  and  3  of  Test  1  occur  in  pairs  and 
involve  two  cells  to  each  side  of  the  cell  boundary  under  consideration. 
From  Figure  9,  Test  11  on  the  above  boundary  includes  this  flux  (GFA) 
and  cells  AA  and  A,  Test  12  includes  GFA  and  cells  B  and  C.  If  one  of 
the  Test  1  conditions  is  true  for  a  boundary,  the  corrector  function 
provides  a  zero  value  of  the  corrected  diffused  conserved  quantity  flux 
for  that  boundary.  If  the  result  of  Test  1  is  false,  the  form  of  the 
corrector  function  is  determined  by  the  sign  of  the  boundary  flux 
functi on . 


In 

terms  of  Figure  6,  if  C  = 

(I.J)  and 

LL  s 

(1-4,0),  etc. 

THE 

VALUE  OF  HF*  IS  ZERO  IF  ANY  OF  THE  FOLLOWING  IS  TRUE. 

1. 

GF*  =  0.0,  *  =  A,B,L,R 

=>  HF*  =  0.0 

2. 

TEST  11 

1 GFA - (DRAA-DRA)l  <  0.0 

=>  UFA  =  0.0 

IGFB- (DRA-DRC)I  <  0.0 

=>  HF  B  =  0.0 

IGFL • (DRR-DRC)I  <  0.0 

=>  HFL  =  0.0 

|GFR - ( DRRR-DRR ) I  <  0.0 

=>  HFR  =  0.0 

3. 

TEST  12 

1 GF  A  - (DRC-DRB)I  <  0.0 

=>  UFA  =  0.0 

IGFB  -  (DRB-ORBB)l  <  0.0 

=>  HFB  =  0.0 

IGFL - (DRL-DRLL)l  <  0.0 

=>  HFL  =  0.0 

1 GFR - (DRC-DRL)I  <  0.0 

=  >  HFR  =  0.0 

where  GF*  -  DIFFUSED  CONSERVED  QUANTITY  FLUX 

DR**  -  DIFFUSED  CONSERVED  QUANTITY  DENSITY 

Figure  9.  ADIFF  phase  part  1  of  the  corrector  function  -  K. 
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Figure  10  contains  the  form  of  the  corrector  function  if 
GF*  >  0.0,  and  Figure  11  the  form  for  GF*  <  0.0.  Use  of  the  asterisk 
character,  *,  is  a  shorthand  notation  for  the  appropriate  index  of  the 
cell  designated  when  the  C,L,R,A,B  letters  are  substituted  for  the 
asterisk.  Use  of  a  double  asterisk,  **,  implies  that  the  second  char¬ 
acter  substituted  indicates  the  direction  relative  to  the  first  char¬ 
acter,  in  which  the  appropriate  variable  is  to  be  found.  If  **  is  set 
to  LL  this  implies  left  of  a  position  to  the  left  of  center,  and  if  ** 
is  set  to  LR,  it  implies  a  position  right  of  a  position  to  the  left  of 
center  which  results  in  the  center  position.  From  this,  the  W  and  P 
functions  can  be  readily  interpreted  in  terms  of  their  independent  var¬ 
iables;  i.e.,  WP*  operates  upon  the  four  diffused  conserved  quantity 
densities  adjacent  to  the  cell  *.  In  terms  of  Figure  6,  WPL  = 

WP  (DRLA,DRLB,DRLL ,DRC) .  The  GF*  and  HF*,  fluxes  are  the  indexable 
boundary  flux  positions  and  correspond  to  the  V F**  positions  in  the 
array  of  Figure  6. 

Figure  12  contains  the  functional  forms  for  the  functions  of 
Figures  10  and  11  which  are  as  yet  undefined.  A  presentation  of  the 
invariant  properties  of  the  various  parts  of  the  corrector  function 
will  introduce  some  order  into  what  is  a  rather  large  collection  of 
variables  and  functions.  First  it  must  be  noted  that  the  basic  three 
kinds  of  independent  variable  are: 

•  The  cell  diffused  conserved  quantity  densities  (DR**) 

•  The  cell  volumes  (Volume**) 

•  The  cell  boundary  diffused  flux  functions  (GF**). 

All  other  quantities  are  composed  of  a  set  of  these  variables.  Each 
corrector  function  (Test  1  is  not  included  in  the  subsequent  discussion), 
is  a  function  of  the  center  cell,  the  four  adjacent  cells,  and  the  four 
central  velocity  functions,  i.e.,  C,A,B,L,R,GFA,GFB,GFL,GFR.  It  is  then 
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In  terms  of  Figure  6,  if  C  =  (I,J)  and 

LL  e  (1-4, J),  etc. 

HF*  =  MAX  (0.0,  MIN  (GF*,  PPLUS*,  PMINUS*)) 

where 

HF*  -  CORRECTED  DIFFUSED  CONSERVED  QUANTITY  FLUX 
GF*  -  DIFFUSED  CONSERVED  QUANTITY  FLUX 
PPLUS*  =  PPLUS  (GF*,  VOLUME*,  DR*,  WP*,  PP*) 
PMINUS*  =  PMINUS  (GF*,  VOLUMEC,  DRC ,  WMC,  PMC) 


and 

VOLUME*  -  VOLUME  OF  CELL,  *  IS  L,R,A,B 
DR*  -  DIFFUSED  CONSERVED  QUANTITY  DENSITY 
WP*  =  WP  (DR*A,  DR*B,  DR*L ,  DR*R) 

WM*  =  WM  (DR*A,  DR*B,  DR*L ,  DR*R) 

PP*  =  PP  (GF*A,  GF*B ,  GF*L ,  GF*R) 

PM*  =  PM  (GF*A,  GF*B ,  GF*L ,  GF*R) 

WHERE  *  IS  C,L,R,A,B,  c.f.  FIGURE  12 

(NOTE  THAT  C*  -  *  ic  CA  =  A,  and  AB  =  C) 


Figure  10.  ADIFF  phase  corrector  function 
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In  terms  of  Figure  6,  if  C  =  (I,J)  and 

LL  =  (1-4, J),  etc. 

HF*  =  (S*)  MAX  (0.0,  MIN  ( |GF*| ,  MPLUS*,  MMINUS*)) 


where 

HF*  -  CORRECTED  DIFFUSED  CONSERVED  QUANTITY  FLUX 
GF*  -  DIFFUSED  CONSERVED  QUANTITY  FLUX 
MPLUS*  =  MPLUS  (GF*,  VOLUMEC,  DRC ,  WP*,  PP*) 

MMINUS*  =  MMINUS  (GF*,  VOLUME*,  DR*,  WM*,  PM*) 

S*  =  S(DR* ,  DRC),  AND  IS  POSITION  DEPENDENT 

AND  THE  REMAINING  FUNCTIONS  AND  VARIABLES  ARE  AS  DEFINED  IN  FIGURE  10; 
i.e. , 

VOLUME*,  DR*,  WP*,  WM*,  PP*,  PM* 

WHERE  *  IS  C,L,R,A,B 


Figure  11.  ADIFF  phase  corrector  function  for  GF*  •  0.0 
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In  terms  of  Figure  6,  if  C  e  (I,J)  and 

LL  =  (1-4,0) ,  etc. 

WP*  =  MAX  (DR*A,  DR*B,  DR*L ,  DR*R) 

WM*  =  MIN  (DR*A,  DR*B ,  DR*L ,  DR*R) 

PP*  =  MAX  (0.0,  GF*L )  -  MIN  (0.0,  GF*R) 
+MAX  (0.0,  GF*B)  -  MIN  (0.0,  GF*A) 
PM*  =  MAX  (0.0,  GF*R)  -  MIN  (0.0,  GF*L) 
+MAX  (0.0,  GF*A)  -  MIN  (0.0,  GF*B) 


SA  =  SIGN  (1.,  (DRA-DRC) ) 
SB  =  SIGN  (1.,  (DRC-DRB) ) 
SL  *  SIGN  (I.,  (DRC-DRL) ) 
SR  =  SIGN  (1.,  (DRR-DRC) ) 


PPLUS*  =  (GF*) • (VOLUME*) • (WP*-DR*)/PP* 
MPLUS*  =  (GF*)- (VOLUMEC)  • (DRC-WPC)/PPC 
PMINUS*  =  (GF*)-  (VOLUMEC)  • (DRC-WMC )/PMC 
MMINUS*  =  (GF*)- (VOLUME*)- (WM*-DR* )/PM* 


Figure  12.  ADIFF  phase  functions  used  for 
corrector  function. 
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additionally  only  a  function  of  the  three  independent  variables  in 
the  direction  of  the  flux  it  is  correcting.  That  is  the  corrector 
for  GFR  includes  the  central  quantities  and  the  quantities  RA,  RB,  RR, 
GFRA,  GFRB,  GFRR.  The  exceptions  are  the  cell  quantities  (1+2,  J+2), 
each  of  which  is  involved  in  two  flux  corrections;  the  diffused  con¬ 
served  quantity  DRAR  is  the  same  quantity  as  DRRA.  The  double  notation 
convention  which  produces  this  seeming  irregulari ty,  is  helpful  when 
the  effects  of  the  mesh  boundaries  is  considered.  It  is  seen  that  the 
central  cell  must  be  at  least  two  cells  distant  from  any  boundary  to 
apply  the  corrector  without  regard  to  boundary  conditions.  This  implies 
that  for  a  considerable  number  of  cells,  the  mesh  boundary  effect  must 
be  included  in  ADIFF  computations:  specifically  any  cell  with  a  mesh 
index  of  i  =  1,2, IQ-1,  IQ,  or  J  =  1,2,JQ-1,JQ,  where  IQ  and  JQ  are  the 
index  limits  within  which  FCT  is  applied. 

One  more  property  of  the  various  corrector  function  is  to  be 
noted.  This  is  that  the  directional  sense  of  each  function  is  inde¬ 
pendent  of  the  direction  of  the  cell  (from  the  central  cell)  upon  which 
its  operation  is  centered.  Any  difference  relation  has  the  same  direc¬ 
tion,  e.g.,  (right-most)-(left-most)  and  (above-most)-(below-most)  for 
PM*,  regardless  of  the  reference  cell  denoted  by  the  *.  Another  inter¬ 
pretation  of  this  property  is  that  the  functions  can  be  translated 
along  an  index  direction;  they  do  not  rotate  their  sense  of  operation 
about  the  central  cell.  Unfortunately,  the  mesh  boundaries  impress  a 
reflective  effect  upon  the  operation  of  ADIFF  when  the  center  cell  is 
in  the  boundary  rows  and/or  columns  of  the  mesh. 
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3.7 


FCT  ANTI-DIFFUSION  -  BOUNDARY  CONDITIONS 


The  guiding  principles  for  treatment  of  ADIFF  when  the  subject 
cell  is  proximate  to  a  mesh  boundary  are: 

•  The  flux  of  diffused  conserved  quantities  is  zero 
across  a  mesh  boundary 

•  The  cells  which  are  OLtside  the  mesh,  in  the  structure 
of  Figure  6,  shall  have  no  effect  upon  the  correction 
of  cell  boundaries  which  are  not  mesh  boundaries. 

The  functions  of  Figures  10,  11,  and  12  do  not  explicitly  represent 
boundary  condition  forms  except  by  imposition  of  the  two  principles 
listed  above. 

The  form  of  the  functions  implies  that  by  assignment  of  special, 
chosen  values  to  the  virtual  cell  and  cell  boundary  quantities;  then  the 
functions  may  operate  correctly  when  a  mixture  of  virtual  and  existing 
cells  and  cell  boundaries  is  to  be  referenced.  This  is  almost  the  case; 
however,  Test  1  requires  a  value  of  the  virtual  cells  which  conflicts 
with  the  value  required  by  the  6F*  greater  than  0.0  and  less  than  0.0 
forms  of  the  corrector  function,  nonetheless  this  method  was  incorporated 
in  the  implementation  with  functions  of  the  indexes  used  in  the  coding 
which  remove  these  conflicts. 

For  Test  1  parts  2  and  3,  of  Figure  9,  coefficient  functions 
for  Test  11  and  Test  12  were  developed  which  multiply  the  relations  by 
a  value  of  +1.0  if  both  cells  of  the  relation  are  mesh  cells,  and 
multiply  +0.0  if  one  (or  both)  of  the  cells  is  virtual,  i.e.,  outside 
the  mesh.  This  method  removes  the  conflict  regarding  specified  values 
for  virtual  cells.  For  example,  application  of  Test  12  is  invalid  for 
cell  R  when  the  left  boundary  of  cell  C  is  a  mesh  boundary,  c.f. 

Figure  6;  cell  L  is  virtual,  and  can  contain  a  specified  value  com¬ 
patible  with  the  corrector  functions  of  Figures  10,  11  and  12. 


There  are  sixteen  possible  boundary  flux  functions  (GF**). 

The  first  boundary  principle  requires  that  the  GF**  functions  be  zero 
when  a  mesh  boundary  is  represented.  Again  a  conflict  is  encountered 
in  the  specified  value  of  a  virtual  cell.  To  insure  that  a  zero  value 
of  GFL  occurs  when  cell  L  is  virtual,  the  conserved  quantity  value  of 
cells  L  and  C  must  be  identical.  This  specification  is  not  compatible 
with  the  corrector  function  operation  on  GFR,  which  requires  cell  L  to 
be  identical  to  cell  R.  To  relieve  this  conflict,  the  index  range  of 
the  routine  which  can  compute  all  sixteen  GF**  values  is  restricted  to 
the  computation  of  GF**  values  which  are  not  virtual  nor  mesh  boundary 
fluxes.  This  index  range  restriction  is  a  function  of  the  index  values 
relating  to  cell  C. 

The  above  two  procedures  allow  specification  of  all  values  of 
virtual  cells  and  virtual  and  mesh  cell  boundary  flux  values,  so  that 
the  more  complex  forms  of  the  corrector  function  may  operate  in  the 
given  functional  forms  regardless  of  the  location  of  cell  C  in  the  mesh. 
When  one  of  the  boundaries  of  cell  C  is  a  mesh  boundary,  specification 
of  the  GF*  flux  of  that  boundary  as  zero,  will  cause  Test  1  -  Step  1 
(Figure  9)  to  be  true  and  the  correction  of  that  flux  is  complete. 

This  allows  specification  of  the  virtual  diagonal  cells  in  this  direc¬ 
tion  from  cell  C,  to  those  values  which  are  appropriate  to  correctors 
of  the  adjacent  cell  C  boundaries.  To  illustrate,  if  the  indexes  of 
cell  C  are  (I  =  1;  J  =  4,  <JQ-1),  then  setting  GFL  =  0.0  determines 
that  HFL  =  0.  Note  that  cells  AL  and  BL  are  not  used  in  this  deter¬ 
mination.  If  cell  AL  is  set  equal  to  cell  AR,  cell  BL  is  set  equal 
to  cell  BR,  and  GFAL  =  GFBL  =  0.0;  then  the  corrector  for  fluxes  GFA 
and  GFB  can  operate  with  the  virtual  cells  AL  and  BL,  and  the  zero 
mesh  boundary  fluxes  GFAL  and  GFBL.  Inspection  of  WP,  WM,  PP,  and  PM 
functions  will  show  that  these  specifications  result  in  no  effect  from 
the  virtual  cell.  Note  that  when  cell  C  indexes  are  (I  =  4,  <IQ-1 ; 

J  =  1)  that  cell  BL  is  appropriately  termed  cell  LB  and  is  set  equal  to 
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cell  LA,  to  not  contribue  to  correction  of  GFL .  When  the  cell  C  indexes 
are  (I  =  2,  IQ-1;  J  =  2,  JQ-1 )  then  the  extreme  cells  are  virtual  and 
must  be  set  equal  to  cell  C,  to  insure  no  effect,  and  the  extreme 
boundary  flux  GF**  is  a  mesh  boundary  and  is  set  to  zero.  It  can  also 
be  shown  that  the  diagonal  cells  will  not  conflict  in  their  required 
specified  values  in  the  mesh  corners;  the  diagonal  cell  is  used  in  only 
one  direction  of  flux  correction,  in  these  cases. 

3.8  FCT  ALGORITHM  IN  HULL 

The  foregoing  discussion  contains  concepts  and  methods  developed 
in  the  process  of  FCT  implementation;  some  items  became  apparent  in  retro¬ 
spect.  Thus,  not  all  of  the  items  are  reflected  in  the  coding.  Several 
storage  structures  and  data  movement  procedures  though  quite  inefficient, 
have  been  included  to  allow  flexibility  in  the  code.  This  facilitates 
debugging,  and  inclusion  of  other  than  major  changes  in  the  form  of  the 
algorithm,  which  may  be  changed  for  improved  results  or  simplifications 
which  are  warranted.  The  FCT  code  as  it  is  currently  constituted  is 
debugged  and  operational,  and  its  absolute  effectiveness  can  begin  to 
be  evaluated.  The  execution  time  overhead  cannot  be  evaluated  with 
this  form  owing  to  the  inefficiencies  which  exist. 

Evaluation  of  the  absolute  effectiveness  is  only  one  consider¬ 
ation  for  inclusion  of  this  form  of  FCT.  If  FCT  in  a  fully  efficient 
form  were  to  provide  numerical  results  comparable  to  a  HULL  run  with 
smaller  cell  size  of  the  same  cost;  the  conclusion  would  be  against 
inclusion  of  additional  code  which  did  not  provide  any  improvement  at 
reasonable  cost.  Thus,  FCT  if  it  proves  to  be  effective  must  also  prove 
to  be  cost  effective. 


To  improve  the  efficiency  of  the  coding  of  the  algorithm 
requires  the  elimination  of  extra  data  movement.  This  is  readily  done 
in  the  present  form.  Such  a  process  al se-.'-'esuJ ts  in  the  addition  of 
more  index  directed  code,  to  which  the  current  form  of  the  code  is 
amenable  as  was  designed.  The  negative  aspect  of  this  type  of  coding 
is  that  the  indexing  is  several  levels  deep  and  becomes  quite  obscure. 
Note  that  the  accessing  implied  in  figure  6  is  two-dimensional,  for 
example  to  use  the  same  line  of  code  to  compute  the  velocity  functions, 
the  indexes  (I,J)  must  be  biased  as  (I+IB,  J+JB)  where  IB  and  JB  are 
unequal  and  have  values  t  1  or  0. 

The  version  of  HULL  used  in  this  implementation  is  managed 
with  the  CDC-UPDATE  system,  and  therefore  the  evaluation  and  coding 
efficiency  efforts  may  proceed  independently  of  each  other. 


SECTION  4 


SUMMARY  OF  RESULTS  OF  FCT  IN  HULL 

This  section  presents  a  summary  of  results  for  a  set  of  two- 

dimensional  HULL  calculations  performed  to  demonstrate  that  the  chosen 

FCT  technique  had  been  implemented  in  such  a  fashion  so  as  not  to 
affect  the  operational  integrity  of  HULL.  These  calculations  also 
provide  a  test  case  to  compare  with  if  AFWL  chooses  to  implement  the 
FCT  technique  in  the  AFWL  version  of  HULL.  Furthermore,  the  existence 
of  the  calculations  allows  a  comparison  between  HULL  with  and  without 
FCT  to  be  made. 

The  test  calculation  chosen  was  a  hypothetical  hot  sphere  of 

2  meters  initial  radius.  The  density  of  the  sphere  was  the  same  as  the 

cold  air  outside,  or  1.225  milligrams  per  cubic  centimeter.  The  sphere 
was  initially  at  rest  and  at  an  elevated  specific  internal  energy, 
which  was  2  x  10^  ergs/gm.  Therefore,  the  initial  pressure  within 
the  hot  sphere  was  about  1  MPa.  The  surrounding  cold  air  was  slightly 

Q 

above  2  x  10  ergs/gm.  These  initial  conditions  produced  a  shock  that 
was  about  0.4  MPa  (i.e.,  about  60  psi )  at  4  milliseconds  after  "detona¬ 
tion."  Figure  13  summarizes  the  initial  conditions,  and  Appendix  B 
presents  various  snapshots  of  the  pressure  and  density  profiles  for 
various  times.  Because  of  the  spherical  symmetry  of  the  chosen  case, 
it  is  expected  that  the  results  for  any  given  technique  will  be  essen¬ 
tially  the  same  along  any  radial  line.  Any  sizable  error  in  the 
difference  technique  in  one  direction  should  show  up  as  an  asymmetry. 
Although  agreement  would  not  prove  that  the  coding  has  been  validated, 
i.e.,  it  is  not  sufficient,  it  is  necessary. 
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Figure  13.  Cylindrical  mesh  geometry,  initial  state, 
time  =  3.355  x  10~3  sec. 
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The  results  show  excellent  agreement  along  the  x  and  y  (radial 
and  altitude  respectively)  directions  as  can  be  seen  by  a  purusal  of 
Appendix  B.  An  example  of  a  representative  case  is  shown  in  Figure  14 
showing  the  results  along  the  y-direction  in  the  first  column  of  cells. 
The  results  along  the  x-direction  in  the  first  row  are  essentially  the 
same.  The  results  differ  along  the  line  where  I=J,  i.e.,  along  the 
radial  which  is  45  degrees  in  elevation.  However,  these  differences 
are  small  and  the  largest  difference  is  seen  to  be  at  the  peak.  The 
large  dot,  labeled  -  see  text,  is  the  peak  along  the  line  where  I=J. 

In  order  to  compare  the  I=J  radial  with  other  radials,  the  index  I  has 
been  multiplied  uy  the  square  root  of  2  so  that  the  cell  indices  are 
proportional  to  the  length  along  the  radial.  (Each  calculation  used 
constant  zoning,  therefore,  the  ratio  of  length  to  cell  number  is  a 
constant. ) 

The  primary  conclusions  of  the  comparison  of  the  HULL  results 
with  the  HULL/FCT  results  are:  (1)  the  implementation  did  not  degrade 
the  operational  status  of  HULL,  (2)  the  implementation  was  accomplished 
with  no  sizable  errors  in  any  >ne  direction  (agreement  with  the  HULL 
without  FCT  implies  correctness  as  well),  and  (3)  the  two-dimensional 
results  with  FCT  are  similar  to  the  AFWL  1-D  results,  where  the  FCT 
technique  seemed  to  improve  spatial  resolution  of  shock  fronts  about  a 
factor  of  two. 

Although  only  a  few  results  are  shown  in  this  report  for  the 
smaller  zone  size  calculation  performed  with  HULL  (i.e.,  problem 
number  1.0005,  Ax  =  Ay  =  5  cm),  the  results  suggest  that  although  FCT 
improves  shock  front  and  contact  surface  resolution,  it  may  not  improve 
the  solution  elsewhere  any  better  than  can  be  done  with  the  same  size 
cell  HULL  calculation.  Basically,  the  2-D  HULL/FCT  results  for  the 
10  cm-zoned  calculation  were  close  to  the  2-D  HULL  results  for  the  same 
zoned  calculation.  The  differences  between  the  2-D  HULL  results  for  the 
5  cm-zoned  calculation  and  the  2-D  HULL  results  for  the  10  cm-zoned  cal¬ 
culation  were  generally  larger  than  the  differences  between  the  two 
10  cm.  calculations  (HULL  and  HULL/FCT). 
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One  final  observation  is  in  order.  At  the  time  of  this  writing 
the  AFWL  implementation  of  FCT  was  not  able  to  work  at  much  higher  over¬ 
pressures.  The  reasons  are  not  known.  Possibly  an  error  exists  in  the 
manner  the  equations  were  formulated,  or  conceivably  there  exists  other 
more  basic  problems.  Clearly,  more  work  needs  to  be  done  in  formulating 
an  acceptable  way  of  using  FCT  with  the  HULL  difference  scheme.  It  may 
be  more  rewarding  to  use  the  fully  multidimensional  FCT  approach,  (see 
reference  8),  without  trying  to  tailor  the  FCT  technique  for  use  with 
the  HULL  difference  scheme.  Since  it  has  been  demonstrated  that  FCT  can 
be  added  to  HULL  in  the  form  recommended  by  AFWL,  it  follows  that  it  is 
possible  to  implement  in  the  simpler  form  as  used  at  NRL. 


SECTION  5 


IMPROVING  REZONE  TECHNIQUES 

After  a  review  of  certain  AFWL  calculations,  SAI  chose  to  imple¬ 
ment  a  computational  subgrid  technique  to  obviate  the  need  for  most 
rezones.  The  computational  subgrid  is  essentially  a  region  that  is 
recognized  by  the  code  to  be  hydrodynamical ly  active,  thus  allowing  it 
to  skip  calculations  elsewhere  in  the  mesh.  Because  of  the  complex 
architecture  which  is  partially  a  result  of  the  many  options  available 
in  HULL,  care  is  required  in  the  use  of  this  technique  just  as  with 
most  rezone  routines.  The  computational  subgrid  technique  is  discussed 
in  some  detail  in  Section  5.1. 

The  remainder  of  the  effort  addressing  this  task  quantified 
the  effects  of  using  various  zoning  algorithms.  In  order  to  develop 
rezone  algorithms,  a  quantitative  understanding  is  needed  of  effects 
caused  by  changes  in  zone  size.  A  rezone  usually  attempts  to  reduce 
the  number  of  zones  in  a  given  problem  by  gradually  increasing  zone 
size.  This  is  done  to  reduce  large  errors  from  use  of  mismatch  in 
zone  size.  These  results  are  discussed  in  Section  5.2. 

5.1  COMPUTATIONAL  SUCGRIDS 

Most  rezone  procedures  in  HULL  were  introduced  as  economy 
measures,  wherein  the  majority  of  the  cells  in  the  mesh  are  maintained 
in  the  hydrodynamical ly  active  region  of  a  calculation.  Then  when  the 
active  region  approaches  the  mesh  boundaries,  the  cell  dimensions  are 
increased  in  size,  the  hydrodynamic  variables  are  averaged  geometrically 
into  the  expanded  mesh,  and  the  calculation  continued.  This  method 
attempts  to  minimize  the  cost  and  execution  time  needed  by  avoiding  the 
overhead  of  processing  excessive  numbers  of  ambient  cells. 
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If  the  HULL  computations  can  be  confined  to  a  subgrid  of  the 
mesh  which  contains  the  hydrodynamical ly  active  cells,  then  the  need 
for  frequent  rezones  in  a  calculation  can  be  avoided.  With  this  method 
a  much  larger  mesh  can  be  used  initially  with  little  overhead,  and  when 
the  mesh  is  maximally  active,  one  of  the  existing  specialized  shock  fol¬ 
lowing  rezones  can  be  invoked,  e.g.,  rezones  3  or  4.*  Many  of  the  cal¬ 
culations  performed  at  AFWL  are  amenable  to  this  computational  subgrid 
form  of  processing,  e.g.,  shock  tubes,  axially  centered  bursts,  applied 
left  or  bottom  boundary  condition  problems.  Figure  15  illustrates  the 
form  this  algorithm  takes  in  typical  two  dimensional  calculations.  In 
the  shock  tube  example,  the  shock  enters  across  the  bottom  boundary 
(J  =  1,  1  <  I  <  IMAX-1).  Three  rows  of  the  mesh  are  hydrodynami- 
cally  active,  and  two  rows  are  allowed  for  propagation  in  the  cur¬ 
rent  cycle.  Therefore,  the  computations  are  limited  to  the  first 
five  rows  J  =  1  to  JQ  in  this  cycle;  rows  JQ  +  1  to  JMAX-1  are 
ambient  and  do  not  change.  If  JMAX  =  100  rows,  this  cycle  costs 
nominally  5%  of  the  cost  without  this  computational  limitation. 

The  ground  burst  example  in  FigurelB  illustrates  the  computational 
subgrid  limitations  in  two  dimensions;  only  four  cells  in  four  rows 
are  processed  in  this  cycle. 

5.1.1  Computational _ Subgrid  Algorithm 

The  algorithm  is  constructed  to  include  cell  (I  =  1,  J  =  1) 
in  the  subgrid  in  all  cases.  This  greatly  reduces  the  overhead  and 
intrusion  of  the  algorithm  coding  in  HULL.  Four  new  variables  are 
required;  two  variables  (IQUIET,  JQUIET)  specify  the  computational 
limits  in  the  current  cycle,  and  two  variables  are  updated  during  the 
cycle  to  force  extension  of  the  subgrid  in  the  next  cycle  (IREZQ, 

JREZQ).  Figures  16  through  20  contain  the  change  deck  for  the  algo¬ 
rithm  in  terms  of  the  SAI  HULL  version  managed  with  CDC  UPDATE. 


*These  are  specific  options  available  in  HULL.  Details  can  be  found 
in  reference  9. 
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Figure  17.  Computational  subgrid  coding,  subroutine  HULUN 
initialization  of  IQUIFT,  JOUIET. 
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Figure  18.  Computational  subgrid  coding,  phase  1  boundary 
routines. 
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Figure  19.  Computational  subgrid  coding. 
Phase  2  boundary  routines. 
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Figure  16  contains  the  algorithm  coding  inserts  and  changes 
in  the  MAIN  subprogram  of  program  HULL.  It  is  seen  that  simple  branches 
are  taken  when  a  cell  index  exceeds  one  of  the  subgrid  limits.  A  cycle 
is  allowed  to  complete  the  row  (5)  loop  for  all  rows  of  the  mesh. 
Thereby,  the  particular  method  used  for  mesh  storage  (disk,  extended 
core,  incore)  and  mesh  input/output  has  no  interaction  with  the 
algorithm,  and  extra  coding  is  not  required  to  ensure  that  mesh 
transfer  from  old  mesh  storage  to  new  mesh  storage  is  complete. 

The  only  phase  which  can  be  performed  in  rows  -  for  which  J>0QU1ET, 
is  Phase  H4,  and  in  particular  only  the  particle  routines  are  active. 
The  multi-material  routines  in  Phase  H4  are  executed  only  within  the 
subgrid  limits. 

Lines  470  through  610  of  Figure  16,  contain  the  intermedi¬ 
ate  variable  update  procedure.  This  code  is  extractive,  and  does  not 
change  any  of  the  computations  in  HULL. 

Figure  17  contains  the  algorithm  inserts  into  subroutine 
HULUN  of  Program  HULL.  The  coding  in  lines  600  through  740  is  also 
required  by  the  cell  by  cell  activity  flag  code  and  the  "*KEEPT0" 

SAIL  directive  for  this  code  is  to  be  set  to  be  compatible  for  values 
of  the  option  "ACTIVE"  other  than  one  and  zero.  Lines  760  through 
950  contain  the  IREZQ,  JRE7Q  initializing  code  which  is  performed  at 
each  calculation  restart;  therefore,  the  four  algorithm  variables  are 
local  to  program  HUM.  and  need  not  be  in  the  "7BL0CK".  The  activating 
conditions  shown  are: 

1.  Presence  of  material  other  than  air 

2.  A  non-ambient  pressure  (11(1)) 

3.  A  non-minimum  absolute  velocity. 

These  conditions  are  sufficient  for  most  calculations  although  other 
conditions  may  be  necessary  in  some  circumstances  including  an  absolute 
specification  of  the  subgrid  limits. 


Figures  18  and  19  contain  the  coping  inserts  for  the  bottom 
(J  =  1 )  and  top  (J  =  JMAX)  boundary  routines.  The  primary  change  is 
substitution  of  IQUIET  for  IMAXM1  (IMAX-1)  in  the  "DO-LOOPS".  The 
IMAXM1  limit  of  the  loops  is  required  for  the  cell  based  activity 
flags  and  for  no  computation  limiting  code  at  all. 

Figure  20  contains  the  IQUIET,  JQUIET  variable  reset  coding 
to  be  entered  into  subroutine  REMESH  of  program  HULL.  This  is  the 
only  place  where  the  subgrid  limits  can  be  reduced  in  value.  The  new 
limits  are  placed  at  the  geometrically  equivalent  cell  in  the  rezoned 
mesh.  Some  extension  is  allowed  because  a  rezoned  mesh  usually  con¬ 
tains  larger  cells,  and  a  two  cell  boundary  between  the  hydrodynamical ly 
active  cells  and  the  cells  beyond  the  subgrid  is  to  be  maintained. 

The  coding  in  Figures  16  through  20  is  untested  at  this  date. 

It  has  been  provided  to  AFWL  with  the  understanding  that  the  initial 
testing  will  be  done  in  one  of  the  two-dimensional  calculations  AFWL 
will  be  running.  However,  the  modifications  in  the  figures  should  not, 
even  if  incorrect,  degrade  the  operational  status  of  the  code  due  to 
the  way  they  are  implemented. 

5.1.2  Limitations  of  the  Computationa 1  Subgrid  Algor i thm 

The  algorithm  as  presented  in  this  memorandum  is  incompatible 
with  the  following  HULL  calculation  options: 

DIMEN  =  3  (three  dimensional  mesh) 

CODE  =  2*  (interactive  dust) 

RAD  0  (radiation  diffusion) 

REZ0NE=  7  (continuous  rezone) 

STRESS  ■  0  (elastic-plastic  stress) 

SURF  :■  0*  (thermal  layer  bottom  boundary). 
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’n«-  'H'f  '  w*h  .in  .r.t-cri sk  (*)  may  or  may  not  be  compatible. 

*  ‘ *"•  •  <  '  •  •  -  *  ,-i*  1 1  1  which  are  outside  the  subgrid  have 

. .  »’•»•>'  I'hm  >s  likely  to  be  compatible.  How- 

'  ■.  ■<  in!  by  gravity  and  will  in  general 

.  •-•••  '  n>  nnly  if  dll  such  particles  are 
’  •  ’  •  ■  "  .••*  •»!.)'  ut«jr i ij  at  each  cycle  is  CODE  =  2 

1  t  '  and  ACTIVE  =  1  is  not  recommended 
■  •  1  '  •  .  minitored.  The  SURF  options  gen- 

•"  i  ■  .  •  •  ••  i  ‘  ■  *  •".(']/  into  several  rows  of  cells  at 

,iv,  if  I  QUIET  is  always  equal  to 

l’!-'-1*-'’  •  ,  i’  ■  If'  i  tlways  greater  than  the  uppermost  row 

w»M  1  "i<  ,  it  ',•(*.  »fien  SURF  0  and  ACTIVE  =  1  are  com¬ 

bat  i:  >  >.  :  r  >  genet  a  i  ,  a  fireball  top  will  force  JQUIET  to 
i nil i .  ife  *  r,.w  We’’  above  the  rows  affected  by  the  SURF  routines,  and 

the  njif  i  on  ;u  n  i  .  i  t n’i| ia  t  ible  . 

sA  I  rei  oiniiends  that  inclusion  of  DEEP-SIX  or  calls  to  SINK  be 
ini  hided  for  the  i  use  where  the  option  ACTIVE  is  equal  to  1  and  incom¬ 
patible  options  are  likely  to  be  accidentally  run.  In  this  manner, 
nonsense  will  not  be  generated  from  some  inadvertent  attempt  to  run 
incompatible  options. 

5.2  QUANTIFYING  THE  REQUIREMENTS  FOR  PROPER  REZONES 

A  series  of  one-dimensional  HULL  (Cartesian  geometry)  calcula¬ 
tions  were  performed  to  study  the  effects  caused  by  changes  in  zone 
sizes  and  zoning  configurations  for  use  in  developing  rezone  techniques. 
Each  run  had  a  reflective  right-hand  boundary  and  identical  initial 
conditions.  The  initial  conditions  consisted  of  air  at  a  density  of 

1.225  milligrams  per  cubic  centimeter  and  a  specific  internal  energy 

g 

of  2.08  x  10  ergs/gm,  which  are  approximately  sea  level  conditions  for 
a  gamma  of  1.4.  The  air  was  initially  at  rest.  Mass,  momentum  and 
energy  were  fluxed  into  the  mesh  from  the  inlet  left-hand  boundary. 
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The  boundary  conditions  were  from  LAMB  (reference  10)  and  represented 
the  air  blast  from  a  1  megaton  nuclear  surface  burst  (2  megaton  free- 
field)  at  either  of  two  ground  ranges,  569  meters  or  8.5  kilometers. 

The  corresponding  peak  overpressures  are  about  4  MPa  and  16  KPa. 

The  waveforms  selected  are  meant  to  approximately  bound  the 
overpressure  region  of  interest.  However,  since  there  exists  interest 
in  investigating  interacting  waveforms,  a  reflective  boundary  was  in¬ 
cluded  at  the  right  hand  side  of  the  mesh  which  is  equivalent  to  the 
case  of  two  equal  shocks  running  head-on  into  each  other.  The  overall 
grid  length  (50m)  and  physical  locations  of  observers  were  essentially 
the  same  from  run  to  run  for  each  given  overpressure.  Table  1  provides 
a  summary  of  the  zoning  and  observer  locations  for  each  calculation. 

The  results  from  each  run  were  processed  with  a  computer  pro¬ 
gram  which  plotted  observer  time  histories  of  various  quantities  of 
interest  (such  as  overpressure).  Figure  21  shows  the  observer  loca¬ 
tions  and  cell  boundaries  for  each  calculation  listed  in  Table  1.  Each 
observer  is  located  at  the  center  of  the  cell  closest  to  the  corres¬ 
ponding  observer  location  in  the  uniform  grid  calculation.  The  symbol 
V  within  parentheses  indicates  that  the  calculation  employed  artificial 
viscosity,  otherwise  no  viscosity  was  used  in  HULL.  The  donor  cell 
technique  in  HULL  will  numerically  generate  some  viscous  dissipation. 
Although  the  treatment  of  the  incident  shock  seems  to  be  handled  nicely 
without  artificial  viscosity,  the  reflected  wave  will  exhibit  a  con¬ 
siderable  overshoot. 

Figure  22  extracted  from  reference  6  shows  three  1-D  HULL  cal¬ 
culation  similar  to  those  in  Table  1  except  that  the  t^flective  wall  is 
at  about  625  m.  The  ones  shown  were  run  with  a  uniform  grid,  one 
without  any  artificial  viscosity,  one  with  the  new  AFWL  artificial 
viscosity  (reference  11),  and  finally  one  with  the  usual  choice  of 
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Table  1.  XO  =  56900  cm 


Nominal  Length  of  Mesh  =  5000  cm. 
(V)  -  Denotes  Vise  Option  with  Cl  =  .6 


RUN 

101 

102 

130,131  (V) 

135,136  (V) 

140,141  (V) 

145,146  (V) 

150,151  (V) 


ZONING 
DX  =  100  cm 

Length  of  mesh  =  5000.0 


DX  =  50  cm 

Length  of  mesh  =  5000.0 


OBSERVERS  &  LOCATIONS 
_ (IQ4  cm) 

1,10,20,30,40,50; 
5.695,  5.785,  5.885, 
5.985,  6.085,  6.185 

1  ,20,40,60,80,100 
same  positions  within 
0.25  m 


0X(1 )  ->  DX (11)  =  100  cm 

DX  =  1.11 24*DX  previous  for  zones  12-26 

Length  of  mesh  =  5001.5  cm 

DX ( 1 )  =  494.2136  cm 

DX  =  (1/1 . 1 1 24 )*DX  previous  for  zones  2-15 
DX ( 1 6)  ► DX ( 26 )  =  100  cm 
Length  of  mesh  =  5001.5  cm 

DX  ( 1 )  -►  DX  (11)  =  100  cm 

DX  =  1.0532*DX  previous  for  zones  12-32 

Length  of  mesh  =  4999.6  cm 

DX ( 1  )  =  296.98  cm 

DX  =  (1/1 ,0532)*DX  previous  for  zones  2-21 
DX( 22) —  DX(32)  =  100  cm 
Length  of  mesh  =  4999.6  cm 

DX  =  100  cm 

Length  of  mesh  =  5000.0  cm 


1,10,17,21,24,26; 
5.695,  5.785,  5.879, 
5.974,  6.076,  6.165 

1  ,3,5,9,16,26; 

5.715,  5.804,  5.876, 
5.981,  6.085,  6.185 


1  ,10,18,24,29,32; 
5.695,  5.785,  5.879, 
5.981,  6.093,  6.175 

1  ,4,8,14,22,32; 
5.705,  5.787,  5.879, 
5.986,  6.085,  6.185 


1  ,10,20,30,40,50; 
5.695,  5.785,  5.885, 
5.985,  6.085,  6.185 
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TABLE  1 .  XO  =  850000  cm 
(continued) 


NOMINAL  LENGTH  OF  MESH  =  100000  cm 
(V)  -  DENOTES  VISC  OPTION  WITH  Cl  =  .6 


RUN 


ZONING 


OBSERVERS  &  LOCATIONS 
(10s  cm) 


132,133  (V)  DX( 1 )  -  DX (11)  =  2000  cm 

DX  =  1.11 24*DX  previous  for  zones  12-26 
DX ( 27 )  =  DX(26) 

IMAX  =  27 

Length  of  mesh  =  100029  cm 


1 ,10,17,21 ,24,26 
8.510,  8.690,  8.878, 
9.067,  9.273,  9.451 


137,138  (V) 


DX( 1 )  =  9884.272  cm  1,3,5,9,16,26 

DX  =  (1/1.1124)*  DX  previous  for  zones  2-15  8.549,  8.728,  8.872, 

DX(16)  =  DX( 27 )  =  2000  cm  9.082,  9.290,  9.490 

IMAX  =  27 

Length  of  mesh  =  100029  cm 


142,143  (V)  DX ( 1 )  -v  DX(ll)  =  2000  cm 

DX  =  1.0532*DX  previous  for  zones  12-32 
DX ( 33 )  =  DX(32) 

IMAX  =  33 

Length  of  mesh  =  999923  cm 


1 ,10,18,24,29,32 
8.510,  8.690,  8.879, 
9.081,  9.305,  9.470 


147,146  'V ) 


DX ( 1 )  =  5939.6  cm  1,4,8,14,22,32 

DX  =  (1/1 ,0532)*DX  previous  for  zones  2-21  8.530,  8.695,  8.878, 

OX ( 22 j  ■>  DX(33)  =  2000  cm  9.092,  9.290,  9.490 

IMAX  =  33 

Length  of  mesh  =  999923  cm 


152,153  (V)  DX  =  2000  cm 

IMAX  =  51 

Length  of  mesh  =  100000  cm 


1 ,10,20,30,40,50 
8.510,  8.690,  8.890,9.090, 
9.290,  9.490 
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OBSERVER 


viscosity  in  HULL.  The  effect  induced  by  the  use  of  artificial  vis¬ 
cosity  is  to  increase  the  numerical  diffusion  as  expected  with  the 
results  that  oscillations  due  to  the  second  order  differencing  are 
damped,  however,  resolution  between  physical  peaks  may  be  lost  if  they 
are  sufficiently  close  together. 

Appendix  C  presents  overpressure  history  for  some  of  the  runs 
and  observers  shown  in  Table  1  and  Figure  22.  To  aid  in  interpretation 
of  the  data,  Figures  23A  and  23B  are  included  here  which  show  the  path  of 
the  various  pressure  peaks  for  the  low  overpressure  and  high  overpres¬ 
sure  uniform  grid  runs  with  no  artificial  viscosity.  The  dots  at  the 
bottom  show  the  locations  of  the  observers  for  convenience.  The  follow¬ 
ing  conclusions  were  made  by  comparing  the  appropriate  overpressure 
histories  that  are  presented  in  Appendix  C: 

1)  The  pulses  are  generally  sharper  and  higher  in 
regions  of  finer  zoning,  and  broaden  and  reduced 
in  regions  of  coarse  zoning. 

2)  The  toe  of  the  incident  pulse  arrives  at  a  given 
observer  slightly  earlier  in  runs  which  include 
the  viscosity,  and  is  especially  noticeable  in 
the  runs  made  at  low  overpressure.  However,  the 
arrival  time  of  the  peak,  when  the  usual  approach 
is  taken  for  defining  peak  arrival  time,  does  not 
depend  on  whether  artificial  viscosity  is  used  at 
either  overpressure. 

3)  The  presence  of  artificial  viscosity  lowers  and 
widens  the  reflected  pulse. 

Other  conclusions  can  also  be  drawn  from  careful  comparison  of 
various  runs  in  the  data  base  represented  by  the  few  examples  shown  in 
Appendix  C.  For  example,  consider  the  comparison  of  low  overpressure 
runs  133,  138,  and  153  at  observer  6.  Each  run  included  artificial 
viscosity,  and  the  uniformily  zoned  run  (153)  produced  the  highest  peak 
and  the  fastest  rise  time.  It  would  seem  that  the  peak  was  reduced 
more  by  propagating  from  small  zones  to  larger  zones  than  vice  versa, 
but  the  rise  times  were  almost  equal  for  the  non-uniform  grid  no  matter 
what  the  direction  of  propagation.  Since  the  peak  measured  at  this 
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Figure  23A.  Low  overpressure  peak  pressure  path.  Run  152 


CALCULATION 

TERMINATED 


REFLECTED  PULSE 


INCIDENT  PULSE 


NO  SHOCKS 

HAVE 

ARRIVED 

2 

3 

4 

5 

• 

• 

• 

• 

i 

579 

1 

589 

1 

599 

- r 

609 

Range  (m) 

Figure  23B.  High  overpressure  peak  pressure  path.  Run  150. 


observer  for  run  133  is  contained  in  a  cell  over  four  times  as  large  as 
that  containing  observer  6  in  both  runs  138  and  153,  the  result  is  not 
surprising.  An  indication  of  the  effect  of  the  cell  size  difference  on 
the  peak  can  be  seen  by  comparing  the  incident  wave  peak  at  observer  1 
for  runs  133  and  138  where  the  observer  is  in  the  small  cell  for  run  133, 
instead  of  138.  The  incident  peak  has  been  reduced  by  a  factor  of  0.83. 
This  is  comparable  with  the  reduction  at  observer  6  of  0.88.  Therefore, 
except  for  the  fact  that  the  distribution  of  pressure  is  over  a  larger 
cell  area,  it  is  expected  that  essentially  the  same  dispersive  effect  is 
seen  when  propagating  from  a  region  of  small  zone  to  larger  zone  and  vice 
versa . 

The  zoning  change  leads  to  pulse  dispersion,  which  can  be  thought 
of  as  a  filtering  effect  where  the  largest  zones  filter  out  frequencies 
higher  than  it  is  capable  of  passing.  This  manner  of  explanation  (i.e., 
use  of  linear  theory)  is  not  appropriate  for  the  high  overpressure,  where 
the  non-linear  effects  are  large. 

Now  consider  the  comparison  of  high  overpressure  runs  131,  136, 
and  151  at  observers  1  and  6.  These  also  were  run  with  artificial  vis¬ 
cosity.  The  uniformily  zoned  run  (151)  produced  the  highest  peak  and 
quickest  rise  time.  Although  the  rise  time  for  run  136  at  observer  1 
was  about  four  times  that  for  run  151  at  the  same  observer,  the  pulse  in 
run  136  at  observer  6  has  sharpened  up  again  by  the  time  it  arrived  due 
to  the  non-linearity  of  hydrodynamics.  Essentially  the  dispersed  pulse 
can  shock  up  since  the  characteristics  within  the  pulse  will  tend  to 
cross  (reference  12). 

This  leads  to  a  significant  and  suggestive  observation  about  the 
problem  associated  with  calculating  low  overpressure  waveforms.  The 
majority  of  experience  in  this  country  in  the  field  of  air  blast  calcu¬ 
lations  has  been  gained  at  the  AFWL.  Generally  the  bulk  of  the  AFWL 
experience  has  been  in  either  computing  fireball  development  (rise, 
growth,  etc.)  or  shock  propagation.  In  the  latter  category,  the  low 
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overpressure  regime  is  the  last  to  be  calculated.  Just  suppose  that  in 
the  process  of  learning  how  to  do  calculations  at  low  overpressure,  the 
experience  gained  at  high  overpressures  guided  the  subsequent  effort. 

It  has  been  demonstrated  in  this  effort,  that  shocks  can  resharpen  at 
high  overpressure  if  the  shock  has  experienced  some  dispersion  as  a 
result  of  propagating  through  regions  of  coarse  zoning.  This  will  not 
happen  at  low  overpressure.  Simply  stated,  what  did  work  as  a  reason¬ 
able  zoning  technique  in  high  overpressure,  will  not  necessarily  work  at 
low  overpressure.  The  implication  is  important  to  the  present  issues 
being  raised  in  the  air  blast  community,  and  therefore,  warrants  addi¬ 
tional  and  more  careful  attention  than  received  to  date.  Figures  24 
and  25,  obtained  by  tracing  the  appropriate  curves  found  in  Appendix  C, 
are  provided  to  show  the  result  for  low  and  high  overpressure.  The  solid 
curves  are  for  the  uniformily  zoned  (lm)  grid,  and  the  dashed  lines  are 
for  the  grid  where  the  zoning  was  coarse  near  the  inlet  boundary  (about 
4m),  and  gradually  (IV  change)  became  finer  as  the  reflective  boundary 
was  approached.  At  observer  6  they  were  lm  wide. 

An  analysis  of  1 -D  HULL  calculations  performed  early  in  this 
effort  had  resulted  in  the  following  interesting  observations  for  high 
overpressures:  The  waveforms  at  observer  6  for  the  runs  where  this 
observer  was  contained  in  a  1  m  cell  were  essentially  the  same  whether 
^he  pulse  has  propagated  through  a  region  which  increased  by  5  or  11 
from  2~  m  t0  1  m  cells  or  even  abruptly,  i.e.,  the  entire  factor  of  2. 
Furthermore,  the  waveforms  at  observer  6  for  the  runs  where  observer  6 
is  contained  in  a  — -  m  cell  were  essentially  the  same  whether  the  pulse 
had  propagated  through  a  region  where  the  zone  size  had  been  decreased 
smoothly  (by  either  5  or  11)  or  abruptly.  The  implication  (at  high 
overpressure)  was  that  the  results  are  fairly  insensitive  to  zoning 
variations.  As  was  just  shown  however,  the  same  was  not  true  of  low 
overpressures.  The  waveforms  at  observer  6  are  shown  in  Figure  26  for 
the  two  uniformily  zoned  cases  where  observer  6  is  in  either  the  small 
(•£  m)  or  large  (1  m)  cell.  The  input  conditions  are  the  same  as  those 
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that  were  discussed  and  shown  earlier.  In  the  figure  oscillations  can 
be  seen  which  result  from  using  the  second  order  HULL  diferencing  (in 
the  Lagrangian  phase)  without  any  artificial  viscosity.  The  period  of 
the  oscillation  is  proportional  to  the  zone  size  at  the  observer  loca¬ 
tion.  Furthermore,  the  rise  time  to  the  peak  is  approximately  equal  to 
the  period  of  this  oscillation.  The  rise  time  is  2  ms  for  the  m  cell 
and  about  4  ms  for  the  1  m  case.  The  rise  time  to  the  pressure  of  the 
incident  pressure  (about  500  psi)  is  about  0.6  ms  for  the  -Jp  m  cell,  and 
since  the  shock  velocity  at  the  incident  pressure  is  about  1.9m/ms  the 
rise  to  the  incident  pressure  is  occurring  over  2-3  cells. 
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Figure  26.  Comparison  of  reflected  pressure  for  m  and  1  n  zoning 


SECTION  6 


CONCLUSIONS  AND  RECOMMENDATIONS 

The  principal  conclusions  of  this  effort  relating  to  the  imple¬ 
mentation  of  FCT  in  HULL  are:  FCT  can  be  implemented  in  HULL  in  a 
reasonable  fashion;  as  currently  configured,  FCT  costs  more  than  merely 
halving  HULL  zone  size;  FCT  cost  can  be  reduced  to  a  more  competitive 
level,  but  coding  becomes  very  complex  and  difficult  to  modify;  and 
more  work  apparently  needs  to  be  done  on  the  FCT  technique  utilized 
at  the  AFWL. 

Our  recommendation  relevant  to  including  FCT  in  HULL  as  formu¬ 
lated  at  AFWL  is  that  the  effort  not  be  pursued  further  since  it  is  too 
risky,  benefits  seem  to  be  minimal,  and  considerable  development  is  still 
required.  Instead,  we  recommend  that  FCT  be  implemented  in  HULL  as 
formulated  at  NRL  (reference  8).  HULL  is  designed  to  allow  easy  modi¬ 
fication  to  the  difference  technique  while  retaining  the  extensive  soft¬ 
ware  needed  for  producing  large  calculations  (e.g.,  appropriate  equation 
of  state  data  bases,  atmosphere  models,  sophisticated  plotting  routines, 
program  tape  libraries,  and  general  problem  data  generators. 

The  principal  conclusions  of  this  effort  relevant  to  improving 
rezone  techniques  are:  use  of  a  computational  subgrid  is  desirable 
since  coding  modifications  are  straightforward  and  in  several  cases  it 
will  obviate  the  need  for  rezone;  zone  size  variations  using  no  arti¬ 
ficial  viscosity  do  not  modify  waveforms  nearly  as  much  as  artificial 
viscosity  variations  for  zone  variations  of  5-10  or  even  abrupt  one¬ 
time  factor-of-2  changes  for  high  overpressure ;  and  effects  from  changes 
in  zone  size  depend  on  shock  strength.  Probably  the  most  important 
observation  has  been  that  pulse  dispersions  introduced  by  propagation 
through  coarse  zoning  will  lead  to  reduction  in  peak  overpressure  and 
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may  never  recover  to  the  correct  value  for  low  overpressures.  Whereas 
for  high  overpressures ,  non-linear  effects  will  tend  to  steepen  any 
dispersed  pulse.  Whether  the  steepening  thus  obtained  produces  an 
accurate  waveform  remains  to  be  demonstrated. 

Our  recommendations  relevant  to  improving  rezone  techniques 
are:  the  SAI  computational  subgrid  technique  should  be  implemented 
and  used  with  the  AFWL  HULL  code;  continue  evaluation  of  effect  of 
zone  size  variation  on  HULL  one-dimensional  results  for  decaying  wave¬ 
forms;  and  finally,  in  order  to  further  refine  rezone  techniques,  the 
effects  introduced  by  use  of  artificial  viscosity  in  HULL  must  be 
better  understood,  as  well  as  the  effects  introduced  by  usage  of  cur¬ 
rent  HULL  rezones. 


> 
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APPENDIX  A 


EXAMPLES  OF  COMPARISONS  OF  CALCULATIONS  WITH  EXPERIMENTAL  DATA 

Reference  4  describes  comparisons  of  theoretical  calculations 
made  with  HULL  and  its  predecessor  SHELL  with  various  high  explosive 
(HE)  experimental  data.  Some  of  the  figures  from  that  reference  have 
been  extracted  and  included  in  this  appendix  for  convenience.  Compar¬ 
isons  of  the  theoretical  calculations  with  high  explosive  (HE)  exper¬ 
imental  data  have  been  performed  at  the  AFWL  as  a  result  of  their 
theoretical  support  of  many  large  scale  instrumented  HE  detonations 
(from  DISTANT  PLAIN  through  DICE  THROW). 

Predictions  of  overpressure  with  SHELL  are  shown  in  Figure  27 
along  with  data  obtained  on  a  500  ton  sphere  of  TNT  detonated  at  the 
ground.  Although  SHELL  results  fall  below  the  experimental  data  below 
pressures  of  10  psi ,  HULL  results  were  considerably  better  due  to  the 
improved  difference  scheme  found  in  HULL.  Experimental  results  are 
plotted  from  PRAIRIE  FLAT  and  DIAL  PACK  which  were  both  detonated  at 
Suffield,  Canada. 

A  comparison  has  also  been  made  for  MINE  UNDER,  which  was  a 
100  ton  TNT  sphere  detonated  1  diameter  above  the  ground.  The  agree¬ 
ment  in  peak  overpressure  between  calculations  and  data  was  similar 
to  that  shown  for  the  tangent  sphere  explosions.  The  positive  phase 
duration  (Figure  28),  which  is  more  difficult  to  measure  than  peak 
overpressure,  and  hence,  has  larger  experimental  error  is  also  shown 
to  be  well  represented  by  calculations. 

Figures  29  through  30  present  data,  taken  by  BRL  and  AFWL 
during  MIXED  COMPANY  (a  500  ton  tangent  sphere  of  TNT),  compared  with 
AFWL  calculations.  Peak  overpressure  agrees  well  as  does  the  over¬ 
pressure  impulse  (Figure  30). 
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MINE  UNDER 


Positive  phase  duration  versus  radius 


COMPANY  (overpressure  versus 


Dynamic  pressure  has  not  been  measured  directly  in  the  past. 

The  normal  method  used  is  to  infer  the  dynamic  pressure  from  measure¬ 
ments  of  the  total  pressure  and  the  overpressure  at  the  same  point. 
Several  smoothing  techniques  are  used  in  an  attempt  to  reduce  any 
accentuated  differences  which  appear  as  high  frequency  oscillations. 

The  inference  of  dynamic  pressure  "data"  for  MIXED  COMPANY  is  further 
complicated  by  the  fact  that  total  pressure  gages  were  not  at  the  same 
locations  as  the  overpressure  gages.  The  former  were  3  feet  above 
ground,  whereas  the  latter  were  at  ground  level.  The  HULL  calcula¬ 
tions  predict  significant  pressure  gradients  between  ground  level  and 
3  feet,  which  if  they  exist  would  tend  to  lower  the  inferred  dynamic 
pressure.  This  may  be  responsible  for  the  disagreement  between  the 
calculations  and  inferred  data  shown  in  Figure  32. 

Reference  4  also  presents  comparisons  with  data  taken  for 
Dipole  West,  shots  8  and  11.  Those  experiments  used  two  charges  deto¬ 
nated  one  above  the  other  such  that  the  distance  between  charges  was 
twice  that  of  the  lower  charge  above  the  ground,  thus  enabling  direct 
comparisons  between  real,  ideal  surfaces,  and  calculations.  Figure  33 
provides  an  example  of  the  excellent  agreement  attained.  From  this 
figure  one  can  infer  that  the  calculations  seem  to  be  providing  the 
flow  field  accurately,  not  just  near  the  ground. 

Some  of  the  most  carefully  performed  experiments  in  recent 
years  were  made  by  Carpenter  at  TRW.  A  large  number  of  triply  redun¬ 
dant  detonations  of  PBX-9404  spheres  were  exploded  over  a  polished 
concrete  slab.  HULL  calculations  were  performed  by  the  AFWL  in  sup¬ 
port  of  this  project.  Figures  34  and  35  show  although  HULL  is  missing 
the  absr  peaks  the  integrated  waveforms  (impulse)  show  excellent 
agreement . 

Finally,  for  PRE  DICE  THROW  Event  2  (a  6  ton  detonation  of  a 
capped  cylinder  of  AN/ FO )  a  HULL  calculation  was  performed  for  a  geom¬ 
etry  never  before  done  and  for  a  relatively  unknown  explosive.  Figure 
36  compares  the  experimental  and  calculated  overpressure  peaks.  Agree¬ 
ment  was  excellent;  the  calculated  waveforms  were  virtually  the  same 
as  those  measured. 


79 


( inches ) 

mind  range  for  8  pound  HE  spheres. 


Ground  range  (inches) 

Figure  35.  Positive  phase  impulse  versus  ground  range  for  8  pound  HE  spheres. 
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APPENDIX  B 

RESULTS  FOR  FCT  IN  HULL 

This  appendix  contains  some  of  the  results  of  calculations  per¬ 
formed  for  the  initial  evaluation  of  the  effects  of  the  Flux  Corrected 
Transport  (FCT)  algorithm  which  has  been  incorporated  into  the  self- 
contained  in-core  version  of  HULL.  This  algorithm  and  its  implementa¬ 
tion  in  HULL  were  described  in  Section  3. 

An  isothermal  sphere  was  selected  as  the  test  calculation,  and 
was  performed  three  times  with  the  conditions  that  follow: 


Calculation  Number 

1 .0010 
1 .0005 
2.0010 


Description 

HULL  without  FCT;  DX=DY^10.  cm 
HULL  without  FCT;  DX=DY=5.  cm 
HULL  with  FCT ;  DX=DY=1 0 .  cm 


In  each  calculation  the  initial  conditions  were  the  same,  and  are  the 
following:  (c.f.  Figure  13). 

Ambient  Conditions: 

-3  -3 

Air  Hass  Density  -  1.22  x  10  gm-cm 

9  -1 

Air  Internal  Energy  -  2.06  x  10  erg-gm 
Isothermal  Sphere: 

Air  Internal  Energy  -  2.00  x  10^  erg-gm 


Sphere  Radius 


-  2.00  x  10  cm. 


During  each  calculation,  the  mesh  variables  in  column  I  -  1, 

Row  J  =  1,  and  diagonal  I  J  J  were  obtained  at  a  set  of  standard  times, 
at  intervals  of  10”^  second.  Figures  37  through  43  contain  histograms 
of  the  pressure  in  the  column  of  cells  adjacent  to  the  Y-axis  (I  -  1,  J) 
c.f.  Figure  13.  Figures  44  through  50  contain  histograms  of  the  mass 


-  •  ,rr; 


•f  ' 


density  in  this  same  column  of  cells.  In  all  histograms  the  solid 
line  represents  data  from  calculation  1.0010  (without  FCT);  the  data 
represented  by  points  are  data  from  calculation  2.0010  (with  FCT). 

The  FCT  algorithm  was  operative  from  the  start  of  calculation  2.0010; 
thus,  these  results  reflect  FCT  throughout  the  duration  of  the  calcu¬ 
lation,  and  do  not  represent  adjustment  of  calculation  1.0010,  at  the 
standard  times  in  the  figures. 

Investigation  of  the  effects  of  FCT  in  one-dimensional  calcu¬ 
lations  were  continuing  at  AFWL.  The  results  presented  in  Figures  37 
through  50  for  the  two-dimensional  isothermal  sphere,  exhibit  behavior 
similar  to  that  in  the  one-dimensional  calculations.  That  is,  some 
clipping  of  the  shock  peak  (of  Figures  37  through  43),  and  a  plateau 
effect  are  observed  in  both  kinds  of  calculations  (reference  6). 

In  its  current  form,  the  FCT  coding  results  in  an  overhead 
factor  of  twenty-five  times  the  cost  of  the  same  calculation  without 
FCT.  This  form  of  the  FCT  coding  is  however  readily  modified.  Some 
work  on  reducing  this  overhead  was  conducted,  and  an  efficient  form 
of  the  FCT  anti -diffusion  phase  can  be  constructed.  The  drawback  to 
implementing  a  more  efficient  form  is  that  the  indexinq  is  many  levels 
deep  and  the  code  is  in  a  form  such  that  almost  any  change  requires  a 
significant  rewrite  and  may  require  another  storage  structure  and 
access  method.  It  was  concluded  that  until  the  final  form  of  the 
algorithm  is  decided,  the  current  FCT  coding  is  adequate  for  FCT  eval¬ 
uation. 

Pressure  and  density  results  along  the  x-axis  were  compared 
with  those  along  the  y-axis.  There  were  no  significant  differences 
between  the  results  along  the  two  axes,  and  therefore  those  along  J=1 
are  not  shown. 
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Figures  51  through  53  contain  plots  of  pressure  at  three 
selected  problem  times.  These  plots  extend  through  the  mesh  along  the 
radial  line  defined  by  the  cells  with  cell  indexes  I  =  J.  The  dis¬ 
tances  between  cell  centers  in  this  direction  is  -J~2  times  greater  than 
that  of  either  the  I  =  1  (column)  or  J  =  1  (row)  directions,  and  there¬ 
fore  these  figures  have  been  plotted  to  the  same  scale  as  the  row  and 
column  plots  by  multiplying  the  cell  index,  I  or  J,  by  J~2.  These  are 
seen  to  be  comparable  to  those  along  1=1. 

Figures  54  through  59  contain  the  pressure  and  density  profiles 
through  the  I  =  1  column  for  calculation  1.0005  at  times  of  4,  5,  and 

_3 

6  x  10  seconds.  The  following  observations  result  from  comparisons 
with  the  corresponding  data  shown  for  the  larger  zoned  calculations . 
Calculation  1.0005  contains  more  detailed  structure  in  the  vicinity  of 
the  shock  front.  The  local  extreme  values  are  sometimes  5  to  10%  dif¬ 
ferent  in  magnitude  when  comparing  the  5  cm  and  10  cm  zone  calculations 
without  FCT.  The  10  cm  zone  calculation  with  FCT  does  not  preferentially 
correct  at  10  cm  zone  calculation  toward  the  values  of  the  5  cm  zone 
calculation. 
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APPENDIX  C 


HULL  I-D  RESULTS 

This  appendix  presents  the  overpressure  histories  for  three 
observers  1,  5,  and  6  as  defined  in  Section  5.  All  runs  shown  here 
were  with  artificial  viscosity.  Two  sets  have  been  included,  the  first 
set  (Figures  59  through  67)  are  for  the  low  overpressure  and  the  second 
(Figures  63  through  76)  are  for  the  high  overpressure.  Figure  21  in 
Section  5  shows  the  zoning  for  each  run.  Individuals  that  desire  the 
waveforms  for  the  other  observers  can  request  same  from  any  of  the 
authors. 
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Figure  61.  Overpressure  run  153  observer 


Figure  63.  Overpressure  run  133  observer  number 


Figure  68.  Overpressure  run  number  138  observer  number 
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Overpressure  run  number  131  observer  number 
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Overpressure  run  number  131  observer  number 
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